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APPLICATION OF LOCAL LINEARIZATION AND 


THE TRANSONIC EQUIVALENCE RULE TO THE FLOW ABOUT 
SLENDER ANALYTICAL BODIES AT MACH NUMBERS NEAR 1.0 
RICHARD W. TYSON 


RALPH J. MURACA 



SUMMARY 


The local linearization method for axisymmetric flow is combined with 
the transonic equivalence rule to calculate pressure distributions on slender 
bodies at free-stream Mach numbers from .8 to 1.2. The method can be used to 
determine the aerodynamic loads on parabolic arc bodies having either 
circular or elliptical cross sections. The method is particularly useful in 
predicting pressure distributions and normal force distributions along the 
body at small angles of attack. The primary application of this method would 
be to establish trends for the flow variables as a function of geometry, and 
to evaluate configurations during the preliminary design stage. The equa- 
tions discussed in this paper may be extended to include wing-body combina- 
tions. 
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INTRODUCTION 


The purpose of this paper is to describe an approximate procedure for 
determining the pressure distributions at free-stream Mach numbers near unity 
for a number of slender bodies both non-lifting and lifting. The quality of 
results will be demonstrated by comparison with experimental data. The 
analysis is based on the small disturbance theory of inviscid transonic flow, 
and makes use of the approximations of slender body theory, the transonic 
equivalence rule, and the method of local linearization. 

Results are presented for several bodies having identical axial distri- 
butions of cross-sectional area and various elliptical cross sections. 

Angles of attack included in these comparisons range from 0° to 6°. The 
examples chosen were selected to enable comparison with existing data 
obtained in wind tunnel tests. 

In addition to a description of the method and a presentation of the 
basic equations used in the derivation of the working equations, the 
computer program listings and input instructions are included as an 
appendix. 
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SYMBOLS 


a 

C 


N 


' N 


M 


dC. 

dx 


N 


Speed of sound 
Pressure coefficient 


P ~ Pco 

Poo 2 

— U 

2 00 


Normal-force coefficient, F^/ qS 


Maximum tody diameter 


Normal force 


m£(y + i) 

__ 

Length of body 

Local Mach number 

Local static pressure 

Free stream static pressure 

Radius of body 

Body cross section area 


Product of reference area and normal force coefficient 
derivative 

Planar oody thickness 
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Free stream velocity 


U 


U l* U 2* 


u. 


u, v, w 
x, r, 9 


0 

a 

Y 

X 

P« 

T 

<t> 



Velocity components 
Perturbation velocity components 

Cylindrical coordinates where x extends in the direction of 
the free-stream velocity 
Total cone angle, in radians 

Angle of attack 

Ratio of npecific heats 

Ratio of major to minor axis length for an ellipse 
Free stream density of air 

Thickness ratio i/d 
Perturbation velocity potential 

Denotes partial derivative with respect to x; 

Denotes second partial derivative with respect 


d± 

9x 

o 

to x; 

3 » 
„ 2 


Subscripts 

00 Free stream conditions 

i, j, k Dummy subscript used in index notation 

0 Stagnation conditions 
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h 


Superscript 


Denotes total derivative with respect to x 
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PRESENTATION OF EQUATIONS 


The equations to be presented here assume very small perturbation 
velocities compared to the free stream velocity U w . The equations of 
motion are then simplified by neglecting small terms in the perturbation 
velocities. The Method of Local Linearization described in ref. 1 is used to 
further simplify the equations and combines with the transonic equivalence 
rule to give approximate solutions to these equations in the transonic flow 
regime. 

Small Perturbation Equations of Transonic Flow 

The equations of motion for steady frictionless flow are given in ref. 2 
as : 


2 5 \ 

i 3Xj a 


3x k 


» J ^ — 1 , 2 , 3. 


( 1 ) 


The speed of sound is related to the velocity components by 



where a Q denotes a reference condition. 

Expanding and substituting the velocity field, defined by 


(2) 


6 



u. » U + u 
1 00 


u 


2 


V 


u 


3 


= w. 


( 3 ) 


where u, v, w are small perturbation velocity components, into equations 
1 & 2 and neglecting terms containing squares of the perturbation velocities 
yields the equation: 


( , m 2 \ 4. 

(1 * V 3^ + 3^ 


= Vt (Y + 


1) — 
U 

oo 


3u 

3x, 


{}- ♦ 

“ 2 


3w 

9x« 


M 2 

M -r ( s; 


3v x 

V 


.,2 w , 3u 

M -iT ( to7 

09 3 


(M 


Retaining terms on the right hand side of eq. ( U ) containing 3u/3x^ , 
and dropping all other higher order terms yields the small perturbation 
equation valid for subsonic, transonic and supersonic flow: 



3u + 3v 
3x 1 3 x 2 


3w . t cL , 
t — = M (y 
3x “ ' 


+ !) H- *3L. 

3x 

oc ]_ 


(?) 


We shall assume the flow to be irrotational and isentropic so that a 
perturbation velocity potential <f> exists, where 


u 



3ij> 

3x 



Thus, the equation for frictionless, irrotational flow with small perturba- 
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tions becomes: 


(1 - »£) 


£$_ , £±_ _ 

< 3 4 *1 ' 


M‘ (Y * 1) 
00 

u~ 


M L . £l_ 

3x l 3x* 


( 6 ) 


Consider a body-fixed Cartesian coordinate system as shown on figure 1, where 
x^ f Xg, correspond to the body axis x, y, z. If the origin is located 
at the nose with the x-axis directed rearward and aligned with the longitu- 
dinal axis of the body, the pressure coefficient can be written (ref 1.); 


C 

P 


U 


* «**> 


where the subscript notation is used to denote derivatives. 


(7) 


The Method of Local Linearization for Axisymmetric Flow 


Approximate solutions of good accuracy for axisymmetric flow with 
1 past a wide class of slender pointed bodies may be obtained by 
application of the method of local linearization. A discussion of the 
development of the method and the resulting equations follow. 

Equation (6) can be rewritten 


where 


X cl) + d> 
xx yy 


<f> = 

zz 



k u 


OO 


( 8 ) 
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and 


k = M 2 (y + 1)/U W 

If X varies sufficiently slowly, it can be considered as a constant in the 
initial stages of the analysis. The equation is then reduced to the 3 - P 
Laplace equation for which a solution is known. This resulting solution l'o,* 
the velocity component u is then improved in the following manner. By 
differentiating the solution for u with respect to x and replacing the 
previously assumed constant X term appearing in the Laplace solution by its 
actual value a non-linear ordinary differential equation results. This 
differential equation may then be integrated numerically. 

This method results in the following first order nonlinear ordinary 
differential equations in u/U w . For M w < 1 the following equation results: 


dtu/uj 

dx 


S ,M (x) 
!*tt 


Jin (l - 


M 2 - 

CO 


k u) + E' (x) 


for 0 < x < Z 
where 


E(x) 


S" (x) 
Lit 


Jin 


S 

inrx(Jl-x) 


+ 



S" (x)-s" U) d€. 
o jx - K\ 


For M m > 1 the following equation results; 


d < u/u J _ 

dx 

where H{ 



Jin (M 2 - 1 + k u) + H' 
00 



(x), 

r s"(x) -s"U) 

J x - E 


dC , 


( 10 ) 


( 11 ) 

( 12 ) 
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for 0< x < i. 

For M w = 1 the following equations result; 
If u < 0 ve have 


d (u/U J _ S'" (x ) 


in (-k u) + P£(x) , 


\ S"(x) „ S(x) .1 rX B" 

p e (x) = TT in ' CTF-x T + FIT -T ““ 


(x) -S“(fJ 

X - f. 


If u > 0 we have 


d (u/U ) 


in (k u) + (x), 


P (x) = S"(x) Jin S(x) 1 f Y S"(x) -S”(C) , r nn 

TT £? *2, o' x-5 « (16 > 

To solve these equations it is necessary to define a value for u 
at some point along the body in order to evaluate the constant of inte- 
gration. The initial value of u may be defined in several different 
ways. If there is no experimental data available for the particular body 
under consideration then there are two classes of cases to be distinguished. 
One class includes the foreparts of bodies with a convex corner, such as a 
cone-cylinder, for which the necessary condition is supplied by the fact 
that the velocity must be sonic at the corner, that is; 

1 - M P 

§ = (IT) 

u * m p (y+D 


s 


The other includes smooth bodies for which a value for u cannot be 
specified a priori. The procedure to be followed in the latter case is to 
find that point (x) along the body which satisfies the condition S"(x) = 0. 
The solution to Laplace's equation at this point will be input as the initial 
u. Numerical integration must progress in both directions from this point. 
The local linearization program which was developed performs this integra- 
tion and provides initial conditions as required. 

“*ie local linearization program solves equations 9 through 16 for 
u/U^ using trapezoidal and th<" Runge-Kutta integration techniques. The 
program then computes pressure coefficients using equation 7. A listing 
of this program and sample input/output are included in the appendix. 

The Transonic Equivalence Rule 
The transonic equivalence rule first stated by Oswatitsch 
(references 3 and U) relates the flow around a slender body of arbitrary 
cross section to that around an "equivalent" nonlifting body of revolution 
having the same longitudinal distribution of cross-sectional area. 

Following the development of reference 1, the general solution 
for the flow over a body can be written 

4 (x, y, z) = (x + az) + 4> (x, y, z) ( 18 ) 

and the pressure coefficient is defined by equation (7). 

From slender body theory: 

4> = 4> 0 + e(x) (19) 

where <|) is the solution of the 2-D Laplace equation for the actual body 
shape, and g (x) is a function of the cross sectional area distribution. 
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Using the equivalence rule the potential can be written 


* = *2,a * *2,t ' *2,B * *B 


(20) 


where is the angle of attacr solution from Laplace's equation, <p~ 

is the thickness solution from Laplace's equation, 4u t, is the solution 

1 13 

to the 2-D Laplace's equation for the equivalent body, <f>_ is the 
equivalent body contribution. 

Consider the application of this method to the case of a body having an 
elliptical cross section. The equivalent body will be a body of revolution 
(circular cross section) having an identical S(x) distribution. 

The equivalent body potential solution is: 


U ® 

^2,B = 2? S ' £n r + 2U ® a r sin 9 » 


( 21 ) 


and 

= - r [<Vb + I? 

The thickness solution is obtained from the potential solution for the case 
of an expanding ellipse. Physically it can be interpreted as the effect 
of passing a body through a plane of fluid normal to the path of motion. 

From reference 5 the solution for a uniformly expanding circle in 
terms of a complex potential, W, is: 

W = r V r log £ (23) 

Where represents the velocity at which the radius of the circle is 

growing, and £ is the complex variable in the transformed plane. The 
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relationship between £ and a , the complex variable in the physical plane 
is 


£ = | [a + (a 2 - a 2 + t 2 ) 2 ] (24) 

where £ » Y + iZ and a and b are the length of the major and minor 
axis of the ellipse. The quantity r Q ^ r can ^e written 


dr dr dx , _ _ 

r at = r to dt = r r u » 


(25) 


but, S’ = 2irr r* 
therefore 


r V - 


S'Ua 

2 7T 


(26) 


The relationship between r, and a and b, is r = -v/afT. Now 
W(£) is the complex potential defined by 


W(0 * ♦ + i* 


where <J> is the perturbation potential in the equation for C . 
Substituting equations 24 and 26 and into equation 23 yields the thickness 
solution for an ellipse 



R.P. 


U S’ 
00 

2 IT 


S-n 


„ . 2 . 2 . 1/2 

g + (q - a + b ) 


(28) 


where R.P. denotes the real part of the potential function. 

Similar analysis of flow over a cylinder (ref. 6) yields the expression 
for angle of attack effects: 
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*2.t * R ‘ P * 


iu «, a r 

— 


(o 2 - a 2 + b 2 


* 1 /2*1 ( a 2 + b 2 ) l/2 

• * 2 2 
J 0 + (o -a + 


0 + (c^-a* 1 + b 2 ) 172 


(29) 


The constants a and b in eq, 28 may be related to the cross sectional 

& 

area of the ellipse, S = irab and a - — . 

Differentiating eq (28) with respect to x yields: 


!k- RP 2= i. 

3x " ‘ ' 2u 3x 


^ + (o 2 - a 2 + b 2 ) 1 ^ 


i0 

Substituting for a, where o = re yieldB, 


%^=K.P. ^ js" tn | ^ 


1 i0 , / i20 , 

le + (e - X + 


+ S’ 


+ An 




i20 , , , n 

e - A + 1/A 


1 \!/2l 
X )\ 

I 


Differentiating with respect to y & z yields, 

<- jt 


94> 2 t U co 

= r.p. — s' 
3y 2 tt 


and 


1 + 


re 


ie 


t 2 i26 2 .2x1/2 

(re - a + b ) 


(30) 


(31) 


(32) 


i9 / 2 i29 2 . 2.1/2 

re + (r e -a +b ) 


3* 2 t U « 

■3-^ = R.P. i— S' 
3z 2 tt 


1 - 


re 


i0 


72 i29 2 ,2.1/2 

(re - a + b ) 


lit 


(33) 


i6 / 2 i29 22.1/2 

re +(r e -a +b ) 



In the angle of attack solution for an ellipse, r 
form the two relationships between a, b and r. 
Differentiating eq 29 with respect to x yields: 








With the derivatives of <p defined all of the terms on the right 

hand side of equation 7 can be evaluated except the pressure coefficient 

distribution on the equivalent body C p . This data must be input from 

E 

either experiment or from the local linearization program previously 
discussed. Evaluation of equation 7 is performed in the Transonic 
Aerodynamics Program. These results may then be integrated to obtain 
values for lift, drag, and pitching moments on slender bodies. Since the 
aerodynamic loading, lift and all lateral forces and moments may be 
expressed in terms of differences in pressure between points at the same 
longitudinal station, these quantities depend solely on A listing of 

the Transonic Aerodynamics Program with sample input/output appears in the 
Appendix . 


APPLICATIONS 

The bodies chosen for analysis fall into three categories, bodies 
of revolution, elliptical bodies, and planar bodies. All were previously 
tested in wind tunnels and experimental data for comparison were readily 
available. Data for four models tested in the Ames lU-Ft Transonic 
Wind Tunnel at Mach numbers from .8 to 1.2 are presented in reference 7. 
Configurations were analyzed and comparisons with experimental data 
are shown in figure 2 through figure 6. 
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The body chosen initially is a body of revolution with axial location 
of maximum cross-sectional area at x/£ = .5, and a fineness ratio of 12. 
The radius of the body is given by 


r(x,8 ) 


Wrr MAX <X - * 2) 

/ p p 2 * 

I cos 9 + A sin 6 


(37) 


where r„. v = 3.0. 

MAX 

Figures 2a and 2b show the pressure coefficients computed by the local 
linearization program, and experimental data from reference 7. The data 
shown on Figure 2a at = 1.0 are in good agreement up to the location 
X./Z = 0.7* where a shock occurs. The results shown on Figure 2b at M » 1.2 
are in good agreement over the entire body. Both cases were run with initial 
conditions obtained from experimental data. Figure 3 shows the differences 
which result when the initial condition is derived from experimental data 
and when it is calculated using the solution at the point where S"(x) = 0. 

The major difference between the results obtained using a computed u^ versus 
an experimentally determined value is a shifting of the curve along the 
axis. 

The second body chosen for analysis was a parabolic arc body of elliptical 
cross-section. The body radius is given by equation 37 above. 
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The Transonic Aerodynamics Program is used to calculate angle of 
attack variation over this body with X - 1. Experimental data was used for 
the equivalent body at zero angle of attack. Comparisons with wind tunnel 
results are shown in figures ka and Vb with Mach number equal to 1.0 and 1.2. 
The pressure coefficient, Cp, is plotted versus x/l at an angle of attack 
of 6°. 

The agreement between computed and measured results is not good for the 
M^ = 1,0 case. The general trend of the data is correct, however, the 
analysis predicts lower value of Cp over most of the body. At - 1.2 
the agreement is better although the analysis continues to underpredict 

c p. 

Figure 5a shows results for a configuration with X = 3 at Mach number 
equal to 1.0. The equivalent body pressure coefficient is again obtained from 
experimental data. The same configuration is shown in figure 5b at 6° 
angle of attack. This analysis was done using the local linearization 
program to predict the equivalent body pressure coefficient used in the 
Transonic Aerodynamic Program. 

For the case at a = 0 the agreement between computed and measured 
data is good. The analysis predicts slightly lower values of Cp over the 
forward portion of the body. Since measured data was used to obtain C p 

E 

the sudden increase in Cp on the rear portion of the body caused by a 
shock is reflected in the computed results. 

For the case at a “ 6° the agreement between computed and measured 
data is not as good. The data at 0 = -90° shows better agreement than that 
at 0 ~ +90°. The general trend in Cp is predicted but the computed values 
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are consistently lower. 

The final body chosen for comparison is the planar body shown in 
figure 6a. The elliptic cone has a major to minor axis ratio of 16.5 and 
is compared to the model tested in reference 1, Shown on figure 6a are the 
computed data for a circular and an elliptical cone with identical cross- 
sectional area distribution. Also shown is experimental data for the 
elliptic cone. As can be seen the agreement is quite good over the forward 
portion of the body; however, as the base of the cone is approached, the 
calculated pressures are higher than measured. This is probably due to wall 
interference effects as cited in reference 1. 

Results for the elliptic cone at a = U° are presented on figure 6b. 
The pressure coefficient distribution along two meridonal locations on the 
body are shown. The data at 0 = +90° represents the windward meridian and 
the 0 = -90° the leeward meridian. In both instances the theory predicts 
higher pressures than measured. The reason for this could be the same as 
cited for the a = 0° results. 


CONCLUSION 

The approximate method presented herein represents a useful tool for 
establishing the variations of pressure and velocity along the surface of 
analytic bodies. The method will appropriately reflect the effects of geometry 
changes on the local pressure coefficients. The computer solutions require 
minimal input and short computation times consequently the method is well 
suited to the analysis of configurations during the preliminary design stages. 
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In this context the quality of agreement between computed and measured results 
is acceptable. The procedures by which these solutions have been obtained 
are not restricted to the particular examples selected for display in this 
paper, but possess much greater generality. 



TION 




□ Measured 



« « 1 . 00 , x = 1.0 

oo 

Figure 3.- Effect of initial conditions on the pressure coefficients 
on a parabolic arc body. 
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Computed 

Measured 



.•I* 


Computed 



I, 


X|o? 




Figure 5(a).- Comparison of pressure coefficients on a parabolic arc 

body with elLiptical cross section at a = 0° and M^- 1.0. 
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APPENDIX 


Local Linearization Program 

Computer Program Input /Output - All input data are input according to the 
NAMELIST type format. In columns 2 through 7 the name $ DATA is punched. 
Column 8 is left "blank. Beginning in column 9 the variable names which are 
shown below and their values are punched. Each variable must be separated 
by a comma and the value of the final variable must be followed by a $ 
sign. If optional variables are not known, they need not be punched. The 
variables may be punched out to column 80 in any order. 

D efinition 

Initial value for X along body 
Pressure coefficient at XO 
/Geometry Variables 
Appearing in body defining 
expression R = A{BX + CX^) 

Ratio oi specific heats 
Mach number 

EXAMPLE: 

$ DATA M = 1.0, XO = .1*8, CP = 0., A = .071, B * 1, 

C = -1, D = 6.03, G = 1.1*, $ 

If more than one case is to be run, then the input data for the second 
case may be placed immediately after the data for the first. 


Optional 


Variable 
XO 
CP 
A 
B | 

C 
D 
G 
M 
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The output consists of a listing of the input parameters, Mach number, 
Initial station and corresponding velocity, and ratio of specific heats for 
the medium. In addition, the equation defining the body shape is shown. 

The variation of velocity v, pressure coefficient CP, and body radius R, 
are also shown versus the longitudinal coordinate X. 


LOCAL LINEARIZATION PROCRAtL 


PROGRAM TAP 1 ( INPUT *0UTPuT#TAPE5= INPUT t PUNCH) 

THIS PROGRAM DETERMINES TRANSONIC AERODYNAMIC CHARACTERISTICS ON 
BODIES of REVOLUTION UR I Nr- THE 
METHOD OF LOCAL LINEARIZATION, 

NAMELIST /OATA/M*X0*U0*A # p,C.D.G,cP 
REAL M*K 

DOOttLE VARC2) »CUVAR(2) 

DIMENSION U<S2)*A(S2)*OER(2> 

COMMON /OER/CUVAR* VARtO£R,M,XOtUO 
COMMON /FC/ A*B*C»DtG»K*Pj 
PI=3. 14159265 

1 XC=100. 

REA04S»OATA) 

1F(E0F*S>2*3 

2 STOP 

3 IF (XO, NE *100, > PRINT 4 
IF (XO.EG.lQO.)PRlNT 5 

4 FORMAT ( lHi i* XO IS INPUT TO THIS PRUORAM, UO IS CALCULATED FROM C 
IP INPUT*) 

5 FORMAT (1H1**X0 AND UO ARE CALCULATED BY THIS PROGRAM*) 
1F(XO.EQ,100.) GO TO 10 

UQ=(-CP-Rl <X0)**2)/2, 

GO TO 11 

CALCULATE XO ANO UO IF NEcESSARY 

10 CALL XUOUO»UO»H) 

11 PRINT 6» M*X0*U0*G»AtB*C»n 

6 FOkMATUH **MACH NUMBER =*tF6,M.5x» "XO =**E16.B,5X» *UQ =**tl6.B/ 

1 IH ♦•GAMMA =*.E16.P// 

2 IH **R (X) = * * FS, 2 , 4h * (*F5.2tilM * <1«X) ♦ ♦F5.2tl0H * (1 
3-X)**,F5.2*lH) ) 

CALCULATE X AND U VALUES ALONG THE BOUy 
CALL XU (U*X*N) 

PRINT OUT X*R*Rl*S*Sl#S?*U*CP 
PRINT T 

7 FORMAT (///lH *l5X»*X**l5X t *U*tl4X,*CP*.15X,*R«.l4X»*Rl*tlAX»*S *♦ 
llAX,*Si*»lAX«*S2*) 

00 B I=1*N 

E=X (1)4 RD*R(E)I R10=Rl(E)4 CPD=CP (U ( I) ♦ E ) 

SDsPl*RlE)**2 

sio- 2 . *pi*h(e)*rHe) 

S20=2,*PI*(R(E)*R2(E) ♦ Ri<E>**2> 

6 PRINT 9 t E#UU) »CPD«R0.R1D,SD*S1D*S20 
9 FORMAT < IH .8EI6.B) 

PLOT Sf S1*S2 VS, X AND U VS, X AND CP VS. X 

CALL PLOT(X.U*N) 

ORIGINAL PAGE IS 
OF POOR QUALITY 


r r c c,c c c 


GO TO 1 
END 

SUBROUTINE XUQ(XQ*U0*M) 

REAL H*K 

01 MENS I ON XI (52) * YI (52) 

COMMON /F C/ A*B*C*0*G*K*Pl 

EPSsI.E-lb 

CALCULATE XO XO=X SUCH THAT S2(X)=0 

XO LIES BETWEEN .02 AND .98 

BISECTION 

X=Q. 

S2X=0. 

1 S2P=S2X 
X=X«.01 

S2X =2.*PI*(R(X )«R2(X ) ♦ RUX >**2) 

1F(ABS(S2X).LT.EPS)G0 TO f- 

IF (S2X*S2P»GE.O.) GO TO 1 

AI=X-.010 

X2=X-.QQ5 

X3 =x 

2 52X1=2. *PI*(R(XI)*r 2(X1) ♦ Rl(Xl)**2) 

S2X2*2.*PI*(R(X2>*R2(X2> ♦ R1(X2)**2) 

IF (ABSCS2XJ .LT*EPS)GO TO (* 

IF (S2 X2 .LT* 0. ) GO TO 4 
IF (S2X1 .LT.O.)GO TO S 
J xl=X2 

X2= (X3*X2)/2. 

GO T02 

4 IF (52 XI .LT .0 . ) GO TO 3 

5 X3~X2 
X2=(Xl*X2)/2. 

GO TO 2 

6 XU=X2 

1F(XO.GE»02.AND.XO.LE.«98)GO TO 8 
PRINT 7 *X0 

7 FORMAT (1H *23H*«* XO OUT cF RANGE ***/ 

1 1H **X0 = **Elb.B) 

STOP 

calculate UQ 

8 IF(M.GE.1.)G0 TO 13 
00 V I=l»Sl 
XI(1)=(I-1.)*.02 

IFUO.LE.XI (I)-.Ol.OR.XO.r-T.XI (I)*.Qi) GO TO 10 
XI (l)sXO 

Yi (1)3-2. *PI*(R(X0)*R3(X0) ♦ 3.*R1 (X0)*R2(/0) > 

GO TO 9 

10 YI (I)s2.*Pl*(R(Xl (I) )*W2(Xl (I) ) ♦ RKXI <I))**2 )/AhS(X0-XI(I)) 

9 CONTINUE 

CALL TRAP(YI*Xl*N»uO) 

GO TO 16 

13 Xl(l)=Q. 

00 14 1=2*51 

14 XI (I)=XI (1-1) *X0/50. 

00 15 1=1 *50 

15 Yl (I)=2.*PI*(R(XI (I) )*H2(xI (I) > ♦ Rl(XKl) >4*2>/(X0-XI<I)> 

YI (S1)=~2.*PI*(R (A0) # R3(XM ♦ J.*p 1 (AO) *H2 <XQ> ) 

N=51 

CALL TRAP(YI,XI*N*UO> 

16 U0=U0/ (4.*3. 14159265) 


Hfe IOHN 
END 

SUBROUTINE XU(U*X.N) 

EXTERNAL DERSU8.CHSUB 
HEAL M*K 

DOUbLE CUVAR(2) ,VAR(2) 

DIMENSION UER (2) *ELEi (2) *fLE2(2) *x<52) *0(52) 

COMMON /OER/CUVAR*VAR # DER»M.XO*UO 

COMMON /F C/ A»3*C»0*G*K»P f 

ELEi(l)=.Q01 

ELE2 ( 1 ) -*001 

K=M**2 # ( 1 • *G) 

INTEGRATE FORWARD 


1 = 0 
11=0 

VAR ( 1 ) =08LE (XO ) 

VAR (2) =DBLE (00 ) 

C1=.0I562S 

1 1=1*1 

CALL INTI (II* 1 * 0 tCI *.02*c»IERR*VAR*CUVAKtDER«tLEl*ELE2tELT*ERRVAL 
1 tDEHSUB tCHSUB*C> 

X ( I ) =SNGL (VAR ( 1 ) ) 

U( I ) =SNGL( VAR(2) ) 

IF (DABS(VARU)-.9«D) .GT..c001O)GO TO 1 

DO 2 0=1*1 

JJ=0-1 

X (5Q-JJ) *X ( I- JJ) 

2 U(SQ-JJ)*U(I-JJ) 

INTEGRATE BACKWARDS 


1=50-1*2 

11=0 

VAR ( 1 ) =DBLE (XO) 

VAR (2) =DBLE (UO) 

CI=-. 015625 
3 1 = 1-1 

CALL INTI (11*1 *0»CI *.02,0.*IERR*VAKfCUVAR*DEH*ELEl*ELE2»ELT*ERRV 
1AL*OERSUB.CHSUB*0) 

X(l)=SNGL(VAR(l)) 

U ( I ) =SNGL ( VAR ( 2) ) 

1F(UABS(VAR(1)-.020).GT..o001D)GO TO 3 

N=S0 

IF (l.EO.l) GO TO 7 
KK=0 

DO 6 J=I *50 
KK=KK*1 
X (KK) =X ( J) 
b U(KK)=U(J) 


N=KK 

7 CONTINUE 
RETURN 
END 

SUBROUTINE DERSUB 


i-VOK oi j 


m r 


rrv 


REAL M*K 

DOUBLE CUVaR (2) * VaR (2) 

DIMENSION DER(2> *XI(52) *Yj (52) 
COMMON /DER/ CUV AH * VARtDEp *M» XO *U0 
COMMON /FC/ A*B*C»D*G*K*Pi 
N=5i 

X=SNGL (CUVAR ( 1) ) 


rr r c r r. 


U=SnGl<CUVAW<?) ) 

S a PI*R(X)**2 

SI « 2.*PI*R<X)«R1<X) 

S d * 2.*Pl*<R(X>«K2(X) ♦ Ri(X)**2> 

SJ a 2.«PI»(R(X)*R3(X) ♦ l.*Rl (X)«*R2(X)> 

$4s2.*Pl*(R<X)*R4(A> ♦ 4. «R l (X) *R3 ( A) ♦ 3.*R2(X)**2> 
IF <H-1 • 1 2*1*7 
1 IF(U.GF.0.) GO TO 7 


M EQUAL TO 1 AND U LESS THAN ' 


calculate integral for usf. in derivative expression 
DELTAsa/50. 

EXTEND LIMITS OF INTEGRATION FOR U .LT.0 
IFU.GT.XO) GO TO 11 
DELTAsXO/SO. 

GO TO 1J 
11 DELTAsX/50. 

13 CONTINUE 
DO 4 1-1*51 
XI U)=(I-1.)*OElTA 

iFTX.GT.Xim-.Ol.AND.X.LT.XniU.Ol) 00 TO S 
S2l=2.»PI*(R(Xl (I))*R2(XI (I) > ♦ Rl (Xl ( 1 ) ) **2) 

VI (I)sHX-Al <I))*5j - S2 ♦ 52I)/(X-XI U) ) »#2 
GO TO 4 
5 Yl(I)=S4 

Xl ( 1 ) =x 
4 CONTINUE 

* CALL TRAP (VI *XI*N*SUH) 

L CALCULATE OERIVaUVE 

DER(2)=(S3*AL06(-K*U) ♦ S 3 ♦ S3*A( OG (S/ <4 ,<*PI*A* < 1 . -X) ) ) 
1 - S2*(1.-2.*A)/(X*U.-X) ) ♦ SUM)/<4.»P1) 

RETURN 

C 

L M ESS THAN I 
l 


C CALCULATE INTEGRAL FOR US,- IN DERIVATIVE EXPRESSION 
2 DO 3 1=1*51 

Xl (I)=(I-1.)*.02 

IF (X .LE*XI (I)-. 01. OR. A .r.T.Xl (I) ♦.01) GO TO F, 
xi d)=x 
vi <n=S4 


GO TO 3 

6 S2I= 2.«PIMR(XI (I))*R2(Xi <I>) ♦ Pl (XI (I))«*2) 

Vi (i)=(ABS(X-Xl (I) )*S3 - *2 ♦ S2I>/(A-XKI>>’>*2 
O Yi(l) = ( (X-XI ( I ) ) 4 S3 - c2 ♦ S2I)/(x-Xl ( I ) ) ’*»2 

3 CONTINUE 

CALL TRAP(YI*XI*N»SUM) 
c CALCULATE DERIVATIVE 

AbAsi, - M**2 - 
IF (ABA.lT.O.) PRINT 20 *ABa 
20 FORMAT ( 1H **ABA = **Elb.Bi 
IF (ABA.lT.O.) ABA=-ABA 

DER(2)=(S3^ALOG( ABA) ♦ S3 # A|,0G (S/ (4.*PI *X* ( l . -X ) ) ) 

1 - S2*U.-2.*X)/(A*(1.-X) > ♦ SUM)/(4.*PI) 

u 1 - S2*(1.-2.*X)/(X*(1.-X) > - SUM ) / ( 4. *P 1 ) 

Q1=S3*AL0G(ABA) 

02=S3*«L0G(S/(4.*PI*X*(1.-X) ) > 

U3=S2*S1/S 

U‘*=-S2*(1.-2.*X)/(X*( 1 .-X) ) 

0 PRINT 129*X»Q1 *Q2*u3*0H*Sl,M 

129 FORMAT ( 1H .bElb.tt) 

RETURN 


S2*S1/S 


S2*S1/S 


«- m greater than i 

D OH M EUUAL TO 1 AND U GREATER THAN OR EUUAL TO 0 


b 


L 


b 


L 


C 


C 




C 


CALCULATE INTEGRAL FOR usr IN DERIVATIVE EXPRESSION 
7 A I ( 1 ) =*Q • 

DELTAsK/SO. 

DO d 1-2»S1 

0 Xi(l)sXIII-l) ♦ DELTA 
DO 9 I=lt50 

S21* 2#*PI*(R<XI (I) )*RJc<Xl (I) > ♦ Rl (Al < 1 ) ) ««2) 
v YlU) = ((X-Xim>*S3 - S2 ♦ S2I)/<X-A1U>>**2 
Y I (51) — SA 

CALL THAP(Yl»Xl#NtSUM) 

CALCULATE DERIVATIVE 
AbAsH*«2 - I. ♦ K*U 
IF (ABA.LT.O.) PRINT 20»ABi 
lF(AdA,LT.O.) ABA=-ABA 

DER (2) =S3 *ALCO( ABa )/<A.*Pi> ♦ 

1 S3 *ALOO(S /<<*.»»PI*X**?) )/(A.«PI) ♦ 

2 52 *<A«Sl -2.*S ) / ( A , *P 1 *A*S ) ♦ 

3 S3 /<2.*P1> ♦ 

A SUM/ ( 2* # P 1 ) 

RETURN 

END 

SUBROUTINE TRAP ( T » A*N» SUM* 

D I Mt. NS l ON X (N) *Y(N) 

SUM=Q* 

DO 1 I - *. » N 

1 SUM=SUM*(Y(I) ♦Y(l-l) )«(X(I >-XU-l> )/2. 

RETURN 

END 

FUNCTION R (X) 

COMMON /FC/ A*8*C.U*G*K*Pf 
EXPONENT CHECK 

1F(U.NE.0«) 60 TO 1 
R =A*<B*<i.-X> ♦ C) 

GO TO 2 

1 R=A* (B*X*C*X**D) 

1 R-A* (B* ( l »“X) ♦ C*U.-A)*«U) 

2 CONTINUE 
RETURN 
END 

FUMCTION R1(X) 

COMMON /FC/ A*B »C *D »6f K »P i 
EXPONENT CHECK 

IF 1D.NE.0..AND.D*NE.1.) Or. TO l 
Rl=-A»B - A*C*D 


GO TO 2 

1 Rl=-A*B-A«C*D*(l»-X)**lO-l .) 

1 Ri =A* (6*C*D*A* # <0-1 •) ) 

2 CONTINUE 
RETURN 
END 

FUNCTION R2(X) 

COMMON /FC/ AtB»C»D»G*KtPl 
^ EXPONENT CHECK 

IF1D.NE.O..AND.O.NE.I..ANO.D.NE.2.) GO TO I 
r2-A*C*0*(D-1,> 

GO TO 2 

1 H2=A»C*D«UD-l,)Ml.-X)**<n-2.> 

1 R2=A # C*D*(D-1.)*x**(D-2.) 

2 CONTINUE 




/ 

KCIUKN 

ENO 

FUNCTION «3(X) 

COMMON /PC/ 4*B»Ctl)tG*N»Pl 
C EXPONENT Ch£CK 

IF <l>.Nt,0*.AND.O.NE.U.ANO*O.NE.H..ANO.U.Nt.3.) GO TO I 
RJ=-A*C*D* (0-1 . > * (L>-2. > 

GO TO 2 

1 P J=A # C < ‘0* (0-1 .) • (0-2 . > *X*« (0-3. ) 

C 1 R 3 =-A*C* 0 *( 0 -l.) # ( 0 - 2 .)*(l .-X)**( 0 - 3 .) 

2 CONTINUE 
RETURN 
END 

FUNCTION R4(X> 

COMMON /PC/ A*B*C*D*G*K*Pt 

c exponent Check 

IF (U.NE.O..ANO.D.NE.1 • .ANR.0.NE.2. .ANU.D.NE.3. .AN0.D.NE.4.) GO TO 
11 

H<»SA«C*0* (0-1 . ) * (0-2 * ) * <0-3. ) 

GO TO 2 

I R 4 =A*C* 0 *( 0 - 1 .)*( 0 - 2 .»*( 0 - 3 .)*X**< 0 -A.) 

U 1 «4sA*C*0*(0-1.)«(D-2.)*<U-3.)»(1.-X)’ > *(0-A.) 

2 CONTINUE 
RETURN 
END 

FUNCTION CP(U*X) 

COMMON /PC/ A*B*C»l)*G*K.Pf 
CPs-2.*U-Rl (X>*«2 

return 
ENU 

bUbKOUTlNt CHSU6 
Rt T urn 
ENU 

SUBROUTINE PLOT (X»U*N) 

DIMENSION X (N) «ll<N> 

Ri I U«N 
END 

*DATA M=1.0*X0 s .A8.CP=0.*As.07l25.B=l »C=-1*D=6.03»G=1 .4» 
iDATA M=l.2*X0=.S6*CP=0.*A^.07l2S.8=l»C=-l*O=6.03*r,= 1.4» 




XO IS INPUT T r THIS P’noiiiiM. i, j I l,A LC'Jl A If 0 F PPM CP INPUT 
MACH KUMRtF *1. *9 = - • »’ l L "* " -l.832tll2 

gamma = i . 4 ‘'o:“~':.:f+c*. 

RCXI = .07 * < l.O' ” (l-X) *■ -l.V~ * '1-Xj** 6. A 3) 

OUTPUT 


X 

2.CC0C00C03-O 2 
4.CCC00000L r -':2 
6,00000000 — 0 2 


1 • 2CC000G0 r -2 I 
1.4CCO9OCP0- A 1 

1 .60000000? *'1 I 

i.ecccooocF-'' i 

2 .OCOOOOCO‘_-0 1 

2 .2 CCOCC CO" - ' L 
2 • 4 CCOOQ COL -O I 


3. CCCC05C1 — "• I 
3. 20000000" -0 l 

3.4 COCCOCT0- * I 
3 *600000 OOF “01 
3.8COOOOCC:-' 1 

4. CCOOGOOO C -0 l 
4.2CCOOOOO r 1 
4.4CCOOOOO C -Ol 
4.6000 JC-nr-' I 
4. 8 C CO 00 CO F -* I 
5. 0CC(J0000 r -0 I 
5.2Ccooooor-: 1 

5 • 4 CGOOOOOf -9 l 
5.6 C COOGCC - ‘ I 

s.eccoooooc-o 1 

6 . A CC 0009P ~ - 1 1 
6.2CC00000L-O I 
6.4CC00CJCC- ' I 
6.6C000000"-0 I 
6.8CC0CC00L-* 1 

7. CCCOOOGOF-0 I 
7.2Crrooonr- 1 1 

7.4CCOOOOOF-01 
7.6CC00CQC:- • I 
7. 8CC00000 : -0 1 

8. CCC0PG9 l '' ,r l 

8. 2C 0000 00*- )I 
8 • 4 COOOO 00 r -0 l 
8. 6CC0CC I 

8. 8 00000 00 P-0 I 

q.occooooi'F->: 1 

9. 2COOOOOO c :-0 I 
9.4CC0C POOF- 9 l 
9 • 6 CCOOOOO: -0 I 

9 . fiCCOOCOOF -* 1 


IJ 

-3. 020 765300-02 
-2.F.42d U9y*:-F 2 
-2. 1 36 213 3 40-02 


-2.5 17 1 24 J 3- -02 
- 2 . 4800 ' 3780-02 
-2.440 j 323 i--»2 
-2.4 If 1081 IF -0 2 
-2. 3732634 J' - ,'2 
-2. 332 36i 760-02 


-2. 9 7 5 ♦ 75 5 90-02 
-1 .97 It it 75E -" 2 
-1.343 -3 91 9 Ur: - 02 
-1 .6 37 >«54 Lr -C 2 
- I • 495 3972 6” ”02 
-1.26 PU 6c"-: 2 
— 4 .75 807 600i —0 3 


4 . 6 r <44 54 020 - 0 2 
5. 3 74 266 dir-; 2 
6. i 3O36470- -02 
t .4 1 77'. 392F-. 2 
7.727o4l 22t-02 

8. 5 4 7-r 7 72 •+-.-* 2 
9.35 C-J 3o 3 l u -02 
I . '• I ?c42 3 l-r ~ • 1 
1.031 3065 5*: -0 1 
1 . I 348 32 3 3--C 1 
I . l6292C45 c -j l 
I • 14 ' I 43 3 OF — 0 1 


4.5 4 789 2 l IF 2 
-4. I 44 } 78965-02 
-1 .49 t 3 56 3 30-51 


CP 

5.53387436E-02 
5.178-.V 8280-92 
4 .9687 743 00-02 

4.8 14 fit 1 760-0? 
4 . 70400668 £ -0 2 
4 .dP 369709 F -;2 
.5 26 703 7bF-02 
4. 452953761 - A 2 
4 • 3bi5.‘ 66 LF-02 
4 .314425060-02 
4 .2413a I 130-02 
4.161 7272 0F-02 
4 ■ C 69i594 76— a 2 
3 .9594 7165E-02 
3 .d26 2ti J50-O2 
3. 662993620-02 
3.462 15 1290-32 
3.2 15344940-02 
2.912543580-02 
2 .541699100-02 
2 .0 89 3 32?.2t-02 
l .537690 7 7F-02 


-9. 342711 1 80-02 
-1 .0 75 766210-': 1 
-1.226073300-01 
-L .3 84 77 7360-0 l 
-1 .850964360-01 
-1 . 7230323 1F-01 
-l .893 305000-0 1 
-2 ■ 72 4 36 4 30—9 1 

~2 • 2 38 4 5 6 7 4 F —0 l 
- 2 . 38521 029 F— C l 
-2 .49471 262 F -0 1 
-2 .537446730-01 


-2 . 181 645060-01 
-1 .502 39 346F -C 1 
5 .243596680-03 
1 .9 7962 7840-0 l 


F 

1.425C0000E-03 
2.849999740-03 
4. 2749969 4f -03 
;• 6 59582 S9E -03 
7.124933510-03 
d . 5 44 8 "0 360 -0 3 
9 • 97 44942 5t -0 ? 

L. 1 39 e 863 tt -C 2 

1.2E2269310-O2 
1.424 565 WE -02 
1.566728340-02 

1. 7CC 69546F-02 
1.850386150-02 
1.5916951 dE -02 
2. 1324931 30-' 2 
2.272606680-02 
2.411343790- ?2 

2. 549958 5 dC -02 
2 . £ 86 66C3 5F-02 
2. 82 160 73 IE -02 

2 . 0 5 4 29 53 40 - “2 
3.08455O58 E-02 
3.211550030 -02 
3.334 754 830-02 
3.4034625 70-02 
3. 566870950-02 

3. 674C7199E -02 
3.774047280-02 
3. 86565 7C90 -02 
3.947631430-02 

4. C1P56C 760 -02 

4.9 76 38 5 4 00-02 
4. 12C836O70-U2 
4 • 1 406 7 24 9 '£ ~C>2 
4. 158 1 725 30 -02 
4. 147121530-02 

4.1 13 5497 O c -02 
4.C53 27C5 of - A 2 
3.964868420 -02 
3. 8446 6 57 10-“2 
3.639309780-02 
3. 495059530-02 
3.2575715 30-02 
? . 5 7 3 786 140 -02 
2.637932460-02 
2.245514C7E-02 
1 .79129347F --2 
1.2696 7 6 4 8 F -02 
6. 746 9 72 360-0 3 


H • C C GGC 0 CP r - ‘ 2 -2.6 6 37 49700“ V 
l.OCOOOOOOr -0 1 -2.6C53J2 )O0-J2 
— 2 .55-31, 52. 'F - .2 


2.6CCOCOC0"-' I -2.2P 50 26 840 -'2 
2 • 8COQOO 00 r -0 1 -2.2235 133 02 -0 2 
- 2 . 1 55824050-02 


-L .26 31201 ) c -f 3 
-l. 3 3 2611290-03 
3 . t 64 56805- -0 3 
7.442002 146-. 3 
1 . 193 3571 30-02 

1.6 72 385 ? 1L-C 2 
2. 1 3 7o52 7uF - 02 
2.7422 52 46! 
3.33 jJ 0J4 7t -02 

3.9 756 694 50 -6 2 


<3 .6 <'t 8 7 6930 -C 3 

9 . 

-9. 7C4 7c 1 7 50-53 
-1 .795599570-02 
-2 .655507880-02 
-3 .575 136610-92 
-4 .564555450-02 
-5.63157 3360-02 
-o. 781300820-02 
-8. C 1755)170 -0 2 


1 . 646/oj8( -p 1 .463 3 555F-I- l 

8.6 7? 7 69 5 5c “0 2 


ORIGINAL PAGE IS 
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TRANSONIC AERODYNAMICS PROGRAM 


Computer Program Input/Output - All input data are input according to 
NAMELIST type format. In columns 2 through 7 the name and INPUT is punched. 
Column 8 is left blank. Beginning in column 9 the variable names which 
are shown below and their values are punched. Each variable must be 
separated by a comma and the value of the final variable must be followed 
by a $ sign. The variables may be punched out to column 80 in any order. 

Variable Definition 


MACH 


Mach number 


CPB(l) 


ITRIG = 0 
= 1 


The array of pressure coefficients for the 
equivalent body of revolution. These values 
are required from X = 0.05 to X = 0.95 in incre- 
ments of .05. Dimension of the array should 
be 17. 

If X f 1, elliptical cross section. 

If X = 1 and body of revolution to be run 
at a # 0 . 


LAM 


Ratio of major to minor axis of ellipse 
X = a/b. 


ALPHA 
THET (1) 

NTHETA 

EXAMPLE: 


Angle of attack in "adians. 

Values of 0 around body at which pressure 
coefficients are desired in degrees. 

Dimension of THET(i) array. 


$ INPUT MACH = 1., CPB(l) = .15, .1, -05, .02, -.01, 

—. 0 ^, —. 06 , —. 085 , — .09 c » - . 1 , —. 1 , — . 1 , —.l, — .1, —. 01 , 

.05, .01, ITRIG = 0, LAM = 3, ALPf = 1.0h7198, 

THET(l) = -90., -hO. , -0., h0., 90., NTHETA = 5, $ 

If more than one case is to be run, then the input data for the second 

case may be placed immediately after the data for the first. 


2h 


All parameters used in this program are complex and in the form 
A + iB. The real or imaginary parts or both may be zero. This reserves 
space for the values when they are not zero. The built-in functions for 
handling floating point are integer computations which have a complex 
counterpart. 


25 


Output 



The output consists of a listing of the input parameters, X, theta, 
and the real and imaginary parts of Cp. 



TRANSONIC AERODYNAMICS _PROGRA?t_ 


PNOOKAM CCP ( INPO I tuUTPUT t TAHk^= INPUT ♦ T APtb=OUTPUt ) 


Kt AL LAM.HACK.LAMi .LAM2.LiiM3.LAM4 .LAM/.LAM 1 * 


COMPLEX A*tt*C»U*E»K»G*n*MA»bH»CP.I It i 2s » 2 ) .El T di 2b » 'd > • ALPh I » C I ♦ 
lCLPo » FLAM t SLAM * t LAM* TL«M t/*Pi»0 l»CTPS(Ji.U.K*U.U»b*T.V 

DIMENSION CPH (cS> »THtT<2bl .blNl <?b> .CubT (2b> » *. (2) 
tUUi VALENCE (rt(H*CH) 

NAMtL IbT/lNPUT/MACh. CPb » L aM ♦ I f p I G » AlP MA » 1 HL T » 1 1 1 HE T A 

CP DATA Wlltt FuP ' aAC t ttOrY 

iIWiu=l Ytb 
ITrtlG=C I’iO 

C i = CMPL X ( 0 • i 1 . 1 
1 WtAU (bfINPUT) 
ik(tUF.b) 100. .1 

3 WKlTt <t>* I nput » 

AUHHA=ALPMA'* , it 1 74bjj 
ALMMi=CMPLA(0 .. AlPhA) 

PAINT 7 .MAbn.L AM ♦ ALPHA 

7 FUHMAT ( lHi . 1SX »«A* . I IX >* f.-fc.TA*»4X, "CP KfcA|"»4A»*CP IMrtG* * 2CX . "MAC* 
1 Nu = «»Elo,S/lH .MAx.^LAl, = ".f-l*-,b/ln *p 4X*"ALPHA = ^»hlb*8//) 

If < i iKlb.LU. 1 ) PKINI J u 

10 FUPMA T ( 1 H to«>A»«CP UATm rijlh AlPna- 0 k OH ' XACT HOIJY*//) 

LAMl=-LAM*i./LAM 
LAM2=LAM/ (LAM* i . ) 

LAMJ=1./(LAM*1. ) 

LAM*t=4.0*U.«LAM)/( l.*LAM» 

KLAH=CmPLX (LAMm ,0 . ) 

SLAM=CMPL A ( LAM i ♦ 0 . ) 

EVALUATE L«-»ITmEIa And t 1 2T‘it TA 
DU 20 J*1*MHFTA 
TMtT ( J) =T hlT (JM. 0174b->3 
THtlA=TrtEI (J) 

IUTmlT=2.0*Thf I < J) 

EH ( J) =CMPla CCUs tTnETtt) *bi N(THt TA) » 

E i 1 2 ( J ) =CMkl X ( CCP ( TO Tml T ) .bl U FUT-.E M ) 
b i N T (J) =biN(THk T t j) ) 

COb T ( J) =CUs ( THE T ( J ) ) 

20 CCMiNUE 

X VmkILS kKUM .OP IU .ob |iv bftPS (.F .ub 5- 1 = CUUNfFP <- OP X 

X — 0 • 

1=1' 

40 X=A*.0b 
i = I ♦ 1 

IF (a.GT J.^l) C-0 10 i 
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L 
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CALCULATE k »k 1 9hd 

kl=K(PklMt) 
kt=K(0OUbLt PWlMtl 

R=.lZ2r>*X 

kl=.122b 

k<J=0. 

ALKbU=ALl)b<k*k> 

Hb(Ji=Kl*k| 
lkbUl=2.#Rbol 
C«bUl=CrtPLA (PSUl *0. ) 
C(kbLl=CMPLA<TKStJi tO.) 
LAM?*-(i:.«K*LAMi> #< * 2 * U.*«-*LAM 3 I**? 
£LAM=CMkLA (LAM 7 » 0. ) 

LAM^=-k*H*LAM + K<*H/L AM 
TLAM=CMHLA(LAMt,,0.> 


DO bO j3l*NTH e lA 

A=L1T( JMCbOPT (F. A IV( J) ♦Fl-'M) 

b=»b*(LlT (J) ♦CbOKt (E l Tc ( J) ♦bLAM) ) 

Ci=-k«SiNT ( J) 

C2=k*CO:jT ( J) 

C=CmHLA (Ci»C2> 


D=sCbuRT (R*k<»E i 1 d A J> ♦LL^M) 

El- (£•*►<* (LAM 2 +L AM j ) )«'■’*<; 
t=LMHLA(El*u.) 

F = CbvJRT (R*k*El T2 ( J) ♦TL«M> 

(jl = k*CUbT ( J) 

02=K*S1NT(J) 

(j = CMPLX (01 ♦ 02 ) 

rll=K"Pc*kI**Kl 
h=CMPLA (Ml * U • ) 

AA = .b*ALHMi# (-Cl“L/L)-t/ (M-t 1 I (J!*i'l-L 4 t/lu«(i(»tll ( J) *u) > ) 

bb=k*Rl*( l ./ (H*E i I ( J> *t >) ■» (Cl *C/K ) 

CCPo=CMPLA (CPH ( l > » •> . > 


CP = ALPrll *kl« (<«./A-A ) —<? • 'MioLLOO ( k^-i ) -lTk^U i - 2, <*k*k 1 <»CL U(^ ( H ) 
l*CkbUl-2.«Al_pHA<* (AA>bb) -(-.b*ALf , rtf-.D*aLkHlMO/U-(t*D*t»M/(U < MK* 
It IT ( J) ♦D) ) ) *R*n l * ( (F ♦(,) /I <> 1 ./ (k«-K i T ( J) »F ) ) > *«V-< AA*rib) 

T mE T A = TmE f ( J ) 

CP=CP*CCPd 


PkluT a 7»A * Th^T ( j) * # ( I > * a ( i'i 
A 7 FukMAKIH .aFI^.m) 
by COwIlMUc 


C 

IF ( MOlj ( I * y ) ) All 1 1 o « <* o 
t"J PhIni 7»MaCh»LAM» AlHmA 
O 0 TO -4 < / 


-friN^UT 

••i«Crt 

^rB 


LmM 

i I wlo 
rt L PHA 
i -»ET 


h | hE 7 A 
•»r, ND 


0.1t*ui* 

0 • I5t *00 * 0*lE*0O* U.5E-01 

-O.bE-Ul* -U.^E-Ol 

-0.1t*u0* -U.it*00* 

0*0* 0.0* 0.0# U • •' ♦ 0.0* 

0 . JE *(j 1 * 


0* 

0* 10**E*00 » 

-Q.15/079r*0l* -Q.tWtn JE*0U» 

0,0* 0*0* 0.0* 0. ' * 0.0* 


0.0* 0*0* 0.0* 0 • f - * 0.0* 

B * 


-0.1L*00* 

O.ot-Ul* 


-0 • U -01 » 
-0.lt* 00* 

u . 1 1 * 0 0 « 


o.o* o.byai j£*ou. o 

0.0* o.o* o.o* u. 

0.0* 0.0* 0.0* 
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”0. lt*00* 

0 , u ♦ 0.0* o.u* 


• 1 D?07^t *01* 0.0* 

0 * 0.0* 0 « u * 


OUTPUT 


A 


ThtT."i CH Kp. aL 


b.UOOQOOQOt-02 

5 -uoouooooe-oc 

S . UOUUUOO 0F.-Q2 
S.UQQUG0G0t>02 
■D.uoouoooot-Oc: 
1 .UOUUOOOOtl-Oi 

i .uooooooot-oi 

1 .UOOOOOOOt-Oi 
1 .UUUUOOOOt-Ol 

1 .ooouoooot-oi 

1 .bOOUOQQOt-Ol 
1 .dOOOOOOOE-Ol 
l.bOOUOOOOfc-Ol 
1 .iCiOUUOOOE-Oi 

1 . DOOUOOOOt-01 

2. U00UQ00QE-01 
2.UOOUOQOOE-01 
2.U00UUQU0E-0 1 
2.UOOUOOOOt-Oi 
?.uOOUOOOOfc-Ol 
2.bOOOOOQOE-01 
a.buouoooofcxn 
a.boouoooot-oi 
a.boouoooot-oi 

2. bUU00QUQE-01 

3. UUUUU000L-01 
3. uOOUOQOOt-Oi 
3.U00U0QU0L-01 
3 •UOOUUGGGL -Gl 
3* UOOUUOOOt-Ol 
3.SOOOOOCOE-U1 
3.b00u0000ti-0i 
3. dOOuOOGOt-G 1 
3.bU0U0000£-Ul 
3.bU0O0OG0t-0i 


-1.S7079Q00L*0m 

-6.y8i3uouot.-oi 

0. 

6,9dl30000t.-0l 
l.b7079000»_*0. 
-1.57O7 VUOOl*0:> 
- 6 .^eU 0000 t-U l 

0. 

6.S>Ml30000t-0 1 
l.b 7 U 7 '* 0 U 0 L 40 <’ 
-i.b70 7yu00t-*0l 

-6.vt»i JUOUOt-O l 

0. 

b.yai jouooc-o i 

1 ,b70 79000 l ♦ Go 
- l»b7U79QOOc.*U- 
-b.yfaiJijOUOt-O I 

0. 

6,yal Joooo h --oi 
l.b7U?9000c.^G i 
-l ,b70 79t)00t*0'' 
-6.yal3ouoot.-oi 

0. 

6.yal J oOOOl-u i 
I. b7u7y000c.«0'< 
-I.b7u7yoo0t-*0' 
-6.vaiJOuoo^-o i 
0. 

6,yai3ouooti-oi 
1 ,b7U 7s»00ut + 0 » 
-1 .b7(J /y ooot *0'. 
-6.ybl3u000c.-0 1 

0. 

6.yaiJuooot-oi 
1 ,b70 7y000L*0 • 


4.b20b242dt-O2 
6. *i6y I 0 2dy£-u2 
1 .AyO‘t726/t-Ul 
1 . 382?2ddot-Ul 

l ,44lJ^b 3 J -U 1 
l.uUbb49bot_-U2 
2.6l974b3ttt-U2 
y.y77iu7/yt-U2 
y. iyo ?372 /t-.-oa 
y . 3 1 1 ob6 -02 

-2.0l 1 .;bb‘*Jt. 
-1.31071 iU4t.-02 

i*. yyo^/^ayt-. -02 

4.4773 /7bbt-02 
b. lb 7 a 1 7ybt -U2 

-4 . 32 4 (i A 1 bUL-02 
-J.323lb3ybt-u2 
1 , 9 blh 03 J 4 r -02 
1 . 6 dl 63 byiE -02 
2. 332S32bOc-U2 
-b. 1 Jbl b6bje -02 

-b.4l4?/b24E-02 
- 1 . 1 Jb l bbbbf*. -02 
-1 . 198Abb3bt-02 
-4. Ab21 263 7L-U3 
-d. GbO l V6b2E-02 
-7.by9btti3At-u2 
-4,26A93l34t:-u2 
-4.164 7«dUbE-G2 
-J.Attlbd3A<rf.-02 
-9,u689(j449t-u2 
-8. o6b3tt23if.-U2 
-6.446431 A 0 h. - 0 2 
-b. 219] 4d22t -02 
-b.D767b3iUl.-02 


CP 1 M no 


b.6043bl9br-02 
3.483241 UOh -02 

1 ,0by2UllUc -02 
b . 03003b 7 7h -02 
b . Db3 7b2tt2F -t.2 
4.622byl 7uh -:)2 
2.1b02bbttJh -02 

2 . 1 364Ab6bF. -fib 

3 . HbU94 1 9bh -02 

4. bb43/49bh -02 

2. /9blbfj92fc*-02 
9 . 2 72 7 J93ttt* -(.3 
9. Ab30tt2y9f r ~03 
2.U197740ofc-02 
2. 7lU*2blHh*-u2 
1 . 13bbdA3/F-02 
1 . 7bb7b4tti F - m 3 
1 . 7/U9J lbh -0 c 
7.1 4bbu 3 1 Of “03 

1 . U4b4bH4 yF' -i ;2 
3.4 14 1 6^6 Jt- -03 
1 . 1 40b j22uh -02 

2. a 607 jJttAh — ;)2 
A .4b66jbl 3c -03 
4 . 293 J/OttbF -0 j 
1 ,b2292&9ut -02 
1 . a H 1 2 3 0 2 1 f -u2 

3. J 3332213c -02 

1 . AbbciJ/y /h “02 

1 . /012J090F -02 
2.0940 1a02F -i)2 

2. b637uhl IF. - 02 
3 . 46Udub6Ah -02 
2. 3U4y lOOUF -02 

2 . /b784AAbF -<j2 


A 


rut r . 




CK IMao 


4.uOOUUOOOL~01 
4.UOOoOOOOf-Oi 
4.gOOuUOOOt-Oi 
4 • OUOuOQOQt-Q 1 
4 , UQGuOOOQt-Q 1 
4.b0000000t-0i 

4.b0ouoooor-oi 
4 .sOOOOOQOE-O 1 
4.:>QQOOGOOfc-Oi 

4. b0ojoooot->oi 

5. UOOUUOOOE-01 
5 . uOQuuoOQt-Ol 
5.00000000E-01 
5.Uu000GQ0t-0i 
5 • OOOOOOOOfc -Ol 
s.boauoooot-oi 
5.b0CO0000E-0i 
5.300000006-01 
5.500000006-01 
5 . toOQUOQOOE-O 1 
b.OQOOOQOOE-Gl 
b.uOQOOQOOE-Oi 
6.000O0000L-01 

b.uooy&ooot-oi 

to.uoouooooe-oi 

b.bOOOOGQOt-Oi 
to. 300000006-01 
fe.bOOOD0006-Ql 
to.bOOOOQQOt-Ol 
to.bGQOOOOOt-Ql 
7 • uOOGGQQOt-Ql 
7.U00000006-01 
7 . UQ 0000006 -0 1 
7 . UUGOCGQOE-Oi 
7.UU0UOG006-01 


-J .57U7M000t-*0 A 
-6. 981300004-01 

0. 

6,98l3u0o0t-0i 
1.570 79UQ0c*0'. 
-I.b7079000c.*0'' 
-b.yeuoooot-oi 
0 • 

6.98iJUUUOt.-Oj 
1 ,370 7V00UL*0 ' 
-l.b7079000t-*0> 
-6.9ai Juooot-o i 

0. 

6.98l30UOUb-0 1 
I.b707v000t.^0^ 
-1.5707 90 GOt.*Q' 
-a. vai JOOOO l-O l 
0. 

to.VfelJOOOOt-Ol 
1 . 570 790004*0'' 
-1.570? 90 00t-*0'' 
-6.981 300004-0 1 
0. 

to.yeuoooot-oj 
1.570790004*0' 
-1 .370790004*0'' 
-6.98l3uQOOL-0 t 

o. 

6.9«1 300004-0 | 
1.3 7u790004*0 > 
-1 .570790004 *Gr> 

-6.vai Juuuot-oi 

0. 

6,4tiiJ00U0t.-0l 
i , 57u 790004*0 .< 


-1 ,0to9S0 7ibt-ol 
-1 . 0 7 1 7b0 1 at-0 1 
-9.1b2 t :892Ufc.-02 
-8.8to333l04t-02 
-8.283^u2bU_-o2 
-l.043l4OtoU.-ul 

-1.065797206-01 
-9.91 b l28oa c -U2 
-9.S99J 37to 76-02 
-4. 10522053c. -u2 
-1 . 0 7‘H044‘*nr.-U 1 
-1. I 1870J4/t-Ui 
-l.120lS9uHt.-Ul 
-I.092837o3t.-Ol 
-1. 054 4 90 JO t. -01 
-1,0244 2b 03t-Ul 
-1 .080*>l4l3t-01 
-1 . ibl 1 33-l‘*t-Ul 
-1 , 1 33?nb99*- i 
-1 . llO-lSJ*r-Oi 
-9 . to 30 2 1 be J t -to 2 
-1 , 051 f-H32>*t-o 1 
-1 . ltoto^bJ4Hr.-Ul 
-1. lto7H4Ue>U.-01 
- 1 . 1 78 •< ? 0 h /t -u 1 
-v.SJb* 5olbt-u2 
-1 . u31 'vl 399t-0 1 
-1.2231 IbOtoL-Ol 

-1 .2494840 7t -Ul 
-1.259 itol*2t-01 
-9«3to5<*0S73£-02 
-1 ,021'->3o2bt“0l 
-1 ,2621 JUo^n -01 
-1 . J21*0UiL-Ul 
-1 ,3b374to»3L-Ul 


-3.to4J2ol 74fc -i?a 
-3. idbtov c 7ut -02 

3 . 7h549«1 7 r. - • » 2 
2. 7«bbU29uf -u2 
3. 30toto2tobbt -02 

-4 . 1 S4b 1 32 1 1 — <j2 

-3.53771*7300* 

3 .8834 /7 inf- -:,ff 
3 . 49U43 1 3o(- -02 
4.1 7b 7 *3uot -08 
-4.b Ih 7O0**U -02 
-3 . 7Q9i 1 3801* —02 
8. 069504 Ihh -ua 
3,ol0to^M43t-02 
4 . d 1 23 7 7 1 7r -f»2 
-h.81uO.JhO 96 "02 
-3 .08999b83t - 11 2 
3 .890924036 — f. 2 
3 .93849b5bf - ,2 

4 , b442i> v23t “02 
-4,42b72379f -02 
-3 . h 7U2vb23h -02 

3. 3bbStobO 7 6 ~uc 
J.«bb7iJ25Jt — 0 2 
4.37869bb3h -02 
-3.9b09b7bot “02 
-3. U39920b9ff -u2 
2.ato?3H78ul-02 
3 . 384 1 a 7536 ”02 
3. o840 74 736 - 02 
-3. i to«94338f -i>2 
-2.38b7o98ut ”02 
2.1 9to8ol i Jr -02 
3. 08380710 r-fi 2 
3.087 JbbO ul “02 
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A 


TMtTA 


O rffcAL 


CH 1 MA 0 


7.bOGOOOOOt-Ql 
7.bOQOOGOOE-Oi 
7 • bQOUOOOOb-Ol 
7 . bQ OUQQOOE-Q A 

7. bUOOOQOOE-Ol 

rt.ooouooooe-oi 

a. uoooooooE-o* 

b. OOOUOOOOt-Ol 
R.UOUOOQQOE-Oi 
H.UOOOOQQQt-Oi 
H.bOOuOOOOE-Ol 

8. b00u(/000t-0i 
8 • bOQOQOQOL -0 i 
8.D00uu000t-01 
8 .bUOuOOOOE — Oi 


-I.b707v000t*0f» 
-6.y«l Joooot-oi 
0 . 

b.veuoooot-oi 

1 ,b70 fvUQOfc.^O'’ 
-l.b 7 O 7 V 0 OOt*O^ 
- 6 .Vdi 3 oOOOC.-Oi 

0. 

6.V81 JOOOOe-Ol 
l.b7C7*000t*0* 
-1 ,S707v000t *0 ^ 
-b.vei jooool-o i 
0 . 

6 .veUuoooi--Oi 

I .b 707 ^ 000 t -*0 l 


-3.222922Vbt-03 
-1.20S3706VL-02 
-<♦.029 32622h-02 
-b,039737<*0L-U2 
-b.6O8iR‘*23t-02 
b.b90>^G31t.-U2 

k, /ObiSVbbot . -02 

l , b^b 9 l 02 bt -02 
3, 2139^1 U /t-0* 

-8. I 36006^11 -03 
1.0371 1 **26l -0 i 

9.b22< I00<«t -02 
b. 1091 3 JH 8 K -02 
3.99 7742U‘*t-02 
2.843631 Jbt-U2 


~2. Jbbbdb2*K-02 
-1 .bObbl^l'+i -92 
l . j<*y^bo3Ht:-u2 
d , Jb24ub2af- -0 d 
1 . v 7b3b i 7af - -*d 
-b.bOV9b^7/t-oJ 

-3.8J9G223of-'-(j J 
3.202* 7^3U(- -0 3 
i .‘♦Obol82Jh -02 
5. 38b7bbl /t* -03 
i . i b 2bl Hi /K --‘i 2 
9. 90 08(’39bf- -J3 
-b.vb 7 b 4 b<*ih,-<ij 

2.078bHb9JF-!jJ 
-1 .2380021 be -02 



